M2.855 · Modelos avanzados de minería de datos · PEC2
2021-1 · Máster universitario en Ciencia de datos (Data science)
Estudios de Informática, Multimedia y Telecomunicación
A lo largo de esta práctica veremos como aplicar distintas técnicas no supervisadas así como algunas de sus aplicaciones reales:
Para ello vamos a necesitar las siguientes librerías:
import random
import umap
import numpy as np
import pandas as pd
import sklearn
from sklearn import cluster # Algoritmos de clustering.
from sklearn import datasets # Crear datasets.
from sklearn import decomposition # Algoritmos de reduccion de dimensionalidad.
# Visualizacion.
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
Este ejercicio trata de explorar distintas técnicas de agrupamiento ajustándolas a distintos conjuntos de datos.
El objetivo es doble: entender la influencia de los parámetros en su comportamiento, y conocer sus limitaciones en la búsqueda de estructuras de datos.
X_blobs, y_blobs = datasets.make_blobs(n_samples=1000, n_features=2, centers=4, cluster_std=1.6, random_state=42)
X_moons, y_moons = datasets.make_moons(n_samples=1000, noise=.07, random_state=42)
X_circles, y_circles = datasets.make_circles(n_samples=1000, factor=.5, noise=.05, random_state=42)
Cada dataset tiene 2 variables: una variable X que contiene 2 features (columnas) y tantas filas como muestras. Y una variable y que alberga las etiquetas que identifican cada cluster.
A lo largo del ejercicio no se usará la variable y (sólo con el objetivo de visualizar). El objetivo es a través de los distintos modelos de clustering conseguir encontrar las estructuras descritas por las variables y.
# NO ME FUNCIONA EN LOCAL (muere el kernel), SÍ EN GOOGLE COLAB
fig, axis = plt.subplots(2, 3, figsize=(13, 8))
for i, (X, y, ax, name) in enumerate(zip([X_blobs, X_moons, X_circles] * 2,
[None] * 3 + [y_blobs, y_moons, y_circles],
axis.reshape(-1),
['Blob', 'Moons', 'Circles'] * 2)):
ax.set_title('Dataset: {}, '.format(name) + ('lo que analizarás' if i < 3 else 'los grupos a encontrar'))
ax.scatter(X[:,0], X[:,1], s=15, c=y, alpha=.3, cmap='jet')
ax.axis('equal')
plt.tight_layout()
En este apartado se pide probar el algoritmo k-means sobre los tres datasets presentados anteriormente ajustando con los parámetros adecuados y analizar sus resultados.
X, y = X_blobs, y_blobs
Para estimar el número de clusters a detectar por k-means. Una técnica para estimar $k$ es, como se explica en la teoría:
Los criterios anteriores (minimización de distancias intra grupo o maximización de distancias inter grupo) pueden usarse para establecer un valor adecuado para el parámetro k. Valores k para los que ya no se consiguen mejoras significativas en la homogeneidad interna de los segmentos o la heterogeneidad entre segmentos distintos, deberían descartarse.
Lo que popularmente se conocer como regla del codo.
Primero es necesario calcular la suma de los errores cuadráticos (SSE) que consiste en la suma de todos los errores (distancia de cada punto a su centroide asignado) al cuadrado.
$$SSE = \sum_{i=1}^{K} \sum_{x \in C_i} euclidean(x, c_i)^2$$Donde $K$ es el número de clusters a buscar por k-means, $x \in C_i$ son los puntos que pertenecen a i-ésimo cluster, $c_i$ es el centroide del cluster $C_i$ (al que pertenece el punto $x$), y $euclidean$ es la distancia euclídea.
Este procedimiento realizado para cada posible valor $k$, resulta en una función monótona decreciente, donde el eje $x$ representa los distintos valores de $k$, y el eje $y$ el $SSE$. Intuitivamente se podrá observar un significativo descenso del error, que indicará el valor idóneo de $k$.
Se pide realizar la representación gráfica de la regla del codo junto a su interpretación, utilizando la librería matplotlib y la implementación en scikit-learn de k-means.
def sse_kmeans(X, k):
sse = []
for i in range(1, k):
kmeans = cluster.KMeans(n_clusters=i)
kmeans.fit(X)
sse.append(kmeans.inertia_)
return sse
def elbow_kmeans(sse, k):
fig, ax = plt.subplots(figsize=(13, 8))
ax.plot(range(1, k), sse)
plt.title("Estimación del número de clústers")
plt.xlabel("Número de clústers (k)")
plt.ylabel("SSE")
plt.show()
# Cálculo
k = 11
sse = sse_kmeans(X, k)
# Visualación
elbow_kmeans(sse, k)
C:\Users\mario\anaconda3\envs\uoc20211pec1\lib\site-packages\sklearn\cluster\_kmeans.py:881: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=4. warnings.warn(
De la gráfica anterior podemos interpretar lo siguiente:
Cabe destacar, que se podría elegir 3 clústers, pero este valor no sería el óptimo ya que la suma de los errores cuadráticos es mayor.
El SSE (suma de los errores cuadráticos), calcula la distancia de cada punto a su centroide, por lo que busca minimizar dicha distancia, es decir, que los clústers sean lo más compactos posibles. Otra forma de mejorar la elección de k es medir la distancia entre centroides, cuanto más distancia haya entre centroides más diferentes van a ser los clústers.
Por lo tanto, se podría hacer el cálculo de maximizar la distancia entre centroides para también mejorar la elección de k, en este caso porque en la anterior ejecución queda claro que lo mejor son 4 clústers, pero hacer este cálculo podría ayudarnos a si es mejor 3 o 5 clústers por ejemplo.
def clusters_view(X, labels):
fig, ax = plt.subplots(figsize=(13, 8))
plt.title("Clústers detectados")
colores = ["blue", "red", "green", "orange", "black", "grey"]
labels_color = []
for i in labels:
labels_color.append(colores[i])
ax.scatter(X[:,0], X[:,1], c=labels_color, alpha=.3)
# Cálculo
modelo_kmeans_blobs = cluster.KMeans(n_clusters=4, random_state=0)
modelo_kmeans_blobs.fit(X)
# Visualización
clusters_view(X, modelo_kmeans_blobs.labels_)
Como resultado de la anterior ejecución, podemos ver que tenemos 4 grupos tal y como era de esperar, ya que al inicializar kmeans el valor de k era 4.
Mencionar que los grupos obtenidos son correctos, es decir, tenemos 4 colores (4 grupos) y los puntos de un mismo color están más próximos al centro del clúster del mismo color que al centro de los otros grupos.
El resultado al ejecutar kmeans nos generá los siguientes grupos:
X, y = X_moons, y_moons
# Cálculo
k = 11
sse = sse_kmeans(X, k)
# Visualación
elbow_kmeans(sse, k)
C:\Users\mario\anaconda3\envs\uoc20211pec1\lib\site-packages\sklearn\cluster\_kmeans.py:881: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=4. warnings.warn(
Según la interpretación de la gráfica el valor óptimo del número de clústers podría ser 4 o 5, en nuestro caso vamos a elegir 4 ya que éste minimiza el error cuadrático de los residuos y generaliza algo más, por lo que a priori esta es la mejor opción.
Al igual que sucedía en el ejercicio anterior, podríamos mejorar la elección de k si calculásemos también la distancia entre los centroides, de tal forma que fuera la máxima, así se conseguría tener clústers bien diferentes unos de los otros a parte de ser homogéneos.
También se podría mejorar la elección de k probando diferentes modelos con diferentes k, en este caso no se sabe a priori si un modelo con 4 grupos es mejor que 5, por lo que se puede generar diferentes modelos y ver los resultado que proporciona.
# Cálculo
modelo_kmeans_moons = cluster.KMeans(n_clusters=4, random_state=0)
modelo_kmeans_moons.fit(X)
# Visualización
clusters_view(X, modelo_kmeans_moons.labels_)
En este caso, vemos que el modelo generado no es bueno, es decir, no consigue clasificar bien las instancias del dataset. Esto se debe principalmente a que el dataset es muy variado.
Aún así vemos que realmente sí que hay grupos, para ser más exactos hay dos, dos medias lunas, pero kmeans no es capaz de detectar bien los grupos básicamente porque el problema ha cambiado, los centroides de cada clúster no están claros, aunque sabemos que hay dos clústers, en resumen es un problema de densidad no de distancia. Kmeans no es capaz de elegir dos buenos centroides porque presupone que los clústers tienen forma convexa, no funciona bien si los clústers son no convexos, es decir, calcula los clústers a partir de la distancia no de la densidad.
X, y = X_circles, y_circles
# Cálculo
k = 11
sse = sse_kmeans(X, k)
# Visualación
elbow_kmeans(sse, k)
C:\Users\mario\anaconda3\envs\uoc20211pec1\lib\site-packages\sklearn\cluster\_kmeans.py:881: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=4. warnings.warn(
De la gráfica anterior podemos concluir que el mejor número de clústers puede ser 4 o 5, ya que es en estos puntos donde se reduce más el error cuadrático de los residuos y así vez conseguimos generalizar. En nuestro caso, se ha decidido elegir un valor de k igual a 4, ya que conseguimos reducir el error y generalizar más.
Al igual que en los apartados anteriores, podemos mejorar la elección de k a partir de maximizar la distancia entre los centroides de cada grupo, o también ejecutando varios modelos y ver cuál nos produce un resultado mejor.
# Cálculo
modelo_kmeans_circles = cluster.KMeans(n_clusters=4, random_state=0)
modelo_kmeans_circles.fit(X)
# Visualización
clusters_view(X, modelo_kmeans_circles.labels_)
En este caso pasa lo mismo que en el ejericio anterior, kmeans calcula los clústers a partir de la distancia de cada punto al centroide y no por densidad, por lo que no es capaz de distinguir bien cuántos clústers necesita y cuáles son éstos.
X, y = X_blobs, y_blobs
# http://exponentis.es/ejemplo-de-uso-de-dbscan-en-python-para-deteccion-de-outliers
# Método que devuelve la regla del codo para un epsilón adecuado a nuestro problema
def elbow_dbscan(X, min_samples):
estimator = decomposition.PCA(n_components = 2)
X_pca = estimator.fit_transform(X)
dist = sklearn.neighbors.DistanceMetric.get_metric('euclidean')
matsim = dist.pairwise(X_pca)
minPts = min_samples # Fijamos el parámetro minPts
A = sklearn.neighbors.kneighbors_graph(X_pca, minPts, include_self=False)
Ar = A.toarray()
seq = []
for i,s in enumerate(X_pca):
for j in range(len(X_pca)):
if Ar[i][j] != 0:
seq.append(matsim[i][j])
seq.sort()
plt.subplots(figsize=(13, 8))
plt.title("Estimación de epsilon")
plt.xlabel("Puntos")
plt.ylabel("Distancia")
plt.plot(seq)
plt.show()
Se ha dictaminado que el min_samples sea 15 ya que era lo que mejor funcionaba con nuestro modeloo, es decir, identificaba bien los 4 clústers, ya que de lo contrario identificaba solo 3. A partir de min_samples se calcula con la función anterior el epsilon estimado.
# Epsilón adecuado
elbow_dbscan(X, 15)
# Cálculo
modelo_dbscan_blobs = cluster.DBSCAN(eps=1, min_samples=15)
modelo_dbscan_blobs.fit(X)
# Visualización
clusters_view(X, modelo_dbscan_blobs.labels_)
En este caso vemos que sí que detecta bien los clústers al igual que antes, con la diferencia de que ahora los clústers se calculan a partir de la densidad y no de la distancia. Con eps indicamos el radio de vecindad de un punto, y con min_samples definimos la cantidad mínima de vecinos que hay dentro del radio de vecindad para considerar el punto dentro del clúster.
Cabe destacar que ahora no todos los puntos pertenecen a un clúster en sí, es decir, hay puntos grises que no forman parte de ningún clúster, en otras palabras tenemos outliers. Estos puntos se caracterizan porque no cumplen tanto la distancia como el número mínimo de vecinos en dicha distancia, por lo tanto, ahora a diferencia de kmeans podemos ver que no todos los puntos deberían pertenecer a un clúster, ya que no presentan las mismas similitudes.
X, y = X_moons, y_moons
En este caso se ha dictaminado que min_samples sea 5 ya que era lo que mejor funcionaba para nuestro modelo (jugando con diferentes valor de min_samples), y a prtir de min_samples obtenemos una estimación del epsilon a utilizar.
# Epsilón adecuado
elbow_dbscan(X, 5)
# Cálculo
modelo_dbscan_moons = cluster.DBSCAN(eps=0.075, min_samples=5)
modelo_dbscan_moons.fit(X)
# Visualización
clusters_view(X, modelo_dbscan_moons.labels_)
A diferencia de kmeans, DBSCAN al basarse en la densidad y no en la distancia en sí, conseguimos que en este caso, es decir, cuando los clústers no son convexos, identifique bien los grupos. Aún así vemos que se sigue habiendo outliers, hay pequeños puntos que no cumplen con el número de vecinos mínimos dentro del radio establecido.
En resumen, al ir calculado para cada punto si dentro del radio epsilon hay un número de vecinos mínimos, conseguimos que se identifiquen bien los clústers, ya que no se tiene en cuenta la distancia en sí de cada punto al centroide como sucede con kmeans, ya que en ese caso el centroide es difícil de establecer.
X, y = X_circles, y_circles
Se ha utilizado un min_samples de 10 porque tras probar diferentes valores era lo que mejor funcionaba, y partir de ese valor se estima el valor de epsilon.
# Epsilón adecuado
elbow_dbscan(X, 10)
# Cálculo
modelo_dbscan_circles = cluster.DBSCAN(eps=0.1, min_samples=10)
modelo_dbscan_circles.fit(X)
# Visualización
clusters_view(X, modelo_dbscan_circles.labels_)
En este caso sucede lo mismo que en el ejercicio anterior, DBSCAN sí que es capaz de identificar bien los clústers ya que se basa en la densidad, es decir, no calcula la distancia entre un punto al centroide del clúster, sino que usa el radio (epsilon) y el número de vecinos mínimo (min_samples) para añadir un punto a un clúster.
Cabe destacar que en este caso también tenemos algunos outliers, pero son muy pocos. En resumen, DBSCAN consigue identificar de forma correcta los clústers.
En este apartado se pide visualizar mediante un dendrograma la construcción progresiva de los grupos mediante un algoritmo jerárquico aglomerativo (estrategia bottom-up). Con ello se pretende encontrar un método gráfico para entender el comportamiento del algoritmo y encontrar los clusters deseados en cada dataset.
X, y = X_blobs, y_blobs
import scipy.cluster.hierarchy as sch
from sklearn.cluster import AgglomerativeClustering
def dendrograma(X, method):
plt.figure(figsize=(13, 8))
plt.title("Dendrograma según " + method)
sch.dendrogram(sch.linkage(X, method = method))
plt.show()
# Dendrograma
dendrograma(X, "single")
# Modelo
modelo_aglomerativo_blobs = AgglomerativeClustering(n_clusters=2,
affinity='euclidean',
linkage='single')
modelo_aglomerativo_blobs.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_blobs.labels_)
# Dendrograma
dendrograma(X, "complete")
# Modelo
modelo_aglomerativo_blobs = AgglomerativeClustering(n_clusters=4,
affinity='euclidean',
linkage='complete')
modelo_aglomerativo_blobs.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_blobs.labels_)
# Dendrograma
dendrograma(X, "average")
# Modelo
modelo_aglomerativo_blobs = AgglomerativeClustering(n_clusters=4,
affinity='euclidean',
linkage='average')
modelo_aglomerativo_blobs.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_blobs.labels_)
Para la realización de este ejercicio se ha calculado con tres tipos de enlace:
Enlace simple: usa la distancia mínima entre elementos. Como se puede apreciar tanto en el dendrograma y los clústers que nos proporiona por salida no es el enlace adecuado, según el dendrograma hay 2 clústers y la ejecución del algoritmo vemos que hay dos clústers también, es decir, se produce el problema de cadena, el cual consiste en la unión de grupos solamente por tener elementos muy próximos entre los mismos.
Enlace completo: usa la distancia máxima entre elementos. En este caso tanto en el dendrograma como en al ejecución del algoritmo nos proporciona los cuatros clústers que estamos buscando. Este tipo de enlace consigue eliminar el efecto cadena pero es sensible a los outliers. Aún así este enlace puede ser bueno para este problema.
Enlace medio: es una mezcla entre el enlace simple y completo trantando de mitigar los problemas del efecto cadena y outliers, tanto en el dendrograma como en la ejecución del algiritmo nos indica de forma clara que lo mejor son cuatro clústers, por lo que este tipo de enlace también pueder ser útil para resolver el problema.
Como nota aclaratoria, el número de clústers ideales según el dendrograma viene determinado por la línea vertical más alta sin ser cortada horizontalmente, de tal forma que en ese línea hacemos un corte horizontal a todo el árbol y nos indica el número de clústers ideales para aplicar el algoritmo aglomerativo.
En resumen, tanto el enlace completo como el enlace medio pueden ser usados de forma idónea para resolver este problema.
X, y = X_moons, y_moons
# Dendrograma
dendrograma(X, "single")
# Modelo
modelo_aglomerativo_moons = AgglomerativeClustering(n_clusters=2,
affinity='euclidean',
linkage='single')
modelo_aglomerativo_moons.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_moons.labels_)
# Dendrograma
dendrograma(X, "complete")
# Modelo
modelo_aglomerativo_moons = AgglomerativeClustering(n_clusters=2,
affinity='euclidean',
linkage='complete')
modelo_aglomerativo_moons.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_moons.labels_)
# Dendrograma
dendrograma(X, "average")
# Modelo
modelo_aglomerativo_moons = AgglomerativeClustering(n_clusters=2,
affinity='euclidean',
linkage='average')
modelo_aglomerativo_moons.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_moons.labels_)
Para la realización de este ejercicio se han seguido usando los enlaces simple, completo y medio.
Respecto al enlace simple, vemos que en este caso el dendrograma nos indica de forma clara que dos clústers son lo ideal (porque la linea azul primera es la más grande y si partimos el diagrama horizontalmente por ahí nos quedan dos grupos). Al ejecutar el algoritmo con dos clústers vemos que realiza una perfecta clasificación, es decir, hay dos claros grupos diferenciados y al hacer uso de la distancia mínima entre elementos del grupo se consigue los clústers que estamos buscando, sin tener el efecto de cadena.
En cuanto al enlace completo, el dendrogrma nos indica que el número de clústers ideales pueden ser dos o cuatro, teniendo en cuenta la línea azul más a la izquierda tenemos 2 clústers, pero si nos fijamos en la línea verde más a la izquierda tendríamos 4 clústers. Como sabemos que lo ideal son dos, al ejecutar el algoritmo vemos que no consigue identificar bien los clústers, debido a que calcula la distancia máxima entre elementos.
Respecto al enlace medio, el dendrograma indica que el número ideal es 2 clústers y al ejecutar el algoritmo vemos que tampoco consigue obtener los clústers que buscamos.
En resumen, para este problema el enlace que mejor funciona es el enlace simple, ya que tiene en cuenta la distancia mínima entre elementos.
X, y = X_circles, y_circles
# Dendrograma
dendrograma(X, "single")
# Modelo
modelo_aglomerativo_circles = AgglomerativeClustering(n_clusters=2,
affinity='euclidean',
linkage='single')
modelo_aglomerativo_circles.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_circles.labels_)
# Dendrograma
dendrograma(X, "complete")
# Modelo
modelo_aglomerativo_circles = AgglomerativeClustering(n_clusters=6,
affinity='euclidean',
linkage='complete')
modelo_aglomerativo_circles.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_circles.labels_)
# Dendrograma
dendrograma(X, "average")
# Modelo
modelo_aglomerativo_circles = AgglomerativeClustering(n_clusters=5,
affinity='euclidean',
linkage='average')
modelo_aglomerativo_circles.fit_predict(X)
# Visualización
clusters_view(X, modelo_aglomerativo_circles.labels_)
Como era de esperar, este problema es similar al anterior por lo que el enlace que mejor funciona es el simple, básicamente porque tiene en cuenta la distancia mínima entre los elementos del grupo y como los clústers están bien diferenciados no se produce el efecto cadena. El dendograma de este enlace nos indica que el mejor número de clústers es 2, es decir, los clústers que estamos buscando.
En cuanto al enlace completo y enlace medio no son idóneos para este tipo de problemas, el primer enlace nos indica que el número de clústers tiene que ser 6 (hay que cortar por la penúltima línea vertical azul), y el segundo enlace nos indica que tiene que ser 5 grupos (hay que cortar por la tercera línea azul). Por lo tanto, los dendrogramas de estos enlaces nos están indicando que no son buenos enlaces para nuestro problema, ya que nos indican que lo ideal son 5 y 6 grupos cuando realmente son 2.
Es posible aplicar una amplia variedad de algoritmos para la reducción de dimensionalidad. Para ello se empleará el dataset MNIST compuesto de miles de dígitos manuscritos del 0 al 9. Donde cada imagen se compone de 784 píxeles (imágenes de 28 x 28), por lo que se parte de un número alto de dimensiones.
X, y = sklearn.datasets.fetch_openml('mnist_784', version=1, return_X_y=True, as_frame=False)
Por lo que cada muestra (las 70k filas del dataset) se componen de 784 dimensiones:
X.shape
(70000, 784)
fig, axis = plt.subplots(1, 10, figsize=(12, 6))
for i, ax in enumerate(axis):
ax.imshow(X[y == str(i)][0].reshape(28, 28), cmap='gray_r')
ax.set_title(str(i))
ax.axis('off')
plt.tight_layout()
Si cada algoritmo obtiene resultados distintos a la hora de reducir la dimensionalidad, ¿qué representación es más fiel a la distribución original?
Antes de reducir las 784 dimensiones originales de cada muestra a 2 para poder visualizarlas en 2 dimensiones, es muy útil conocer, o al menos intuir, la estructura en alta dimensionalidad de los datos.
Para ello se puede hacer uso del dendrograma como heurística para conocer la disposición original de los datos y comprobar si la proyección es similar a lo mostrado por el dendrograma.
# Creación del PCA
modelo_pca = decomposition.PCA(n_components=2)
X_pca = modelo_pca.fit_transform(X)
# Visualización
# https://www.kaggle.com/parulpandey/part1-visualizing-kannada-mnist-with-pca
y_numerico = [int(i) for i in y] # Convertimos a número las etiquetas
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(X_pca[:, 0], X_pca[:, 1], s=5, c=y_numerico, cmap='Spectral')
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('Análisis de componentes principales', fontsize=14)
plt.xlabel('Primer componente principal', fontsize=14)
plt.ylabel('Segundo componente principal', fontsize=14)
Text(0, 0.5, 'Segundo componente principal')
Lo primero de todo, las clases no han quedado visiblemente separadas del todo, es decir, sí que se pueden apreciar patrones gracias a los colores, pero están todos los puntos juntos y es difícil de diferenciar en la parte central de la representación.
Uno de los motivos por los que no queda del todo bien diferenciadas las clases puede ser porque estamos reduciendo las dimensiones originales (784) a dos dimensiones, es decir, el modelo no es capaz de representar toda la variabilidad del mismo ya que hay mucha diferencia entre las dimensiones originales a dos dimensiones (se pierde información). Ésto se podría solventar haciendo un estudio de cuántas dimensiones son ideales para este problema, y a partir de ahí generar el modelo con dichas dimensiones. Otro punto pueder por que haya demasiada información en el dataset, es decir, demasiados registros.
Aunque se ha comentado que es difícil de interpretar en profundiad la representación, sí que podemos ver patrones en los datos. Por ejemplo, vemos que el número 0 (rojo oscuro) y el número 7 (verde azulado) están opuestos el uno del otro, esto tiene toda la lógica del mundo y es que los píxeles usados para representar ambos son muy diferentes. Otro caso sería el del 1 (rojo) junto con el 0 y el 7, es decir, el uno usa píxeles muy parecidos al número 7 por lo que ambos están más juntos en la representación que respecto al número 0. También vemos que en el medio de la representación hay muchos números que usan píxeles parecidos, como el 4, el 5, el 6.
En resumen, aunque las clases no quedan del todo definidas sí que lo hacen de forma suficiente como para poder detectar patrones en los datos.
En la gráfica anterior cada punto representa una muestra en 2 dimensiones. Con PCA es posible invertir la transformación para que, a partir de cada punto 2d, se obtenga de nuevo (aproximadamente) la imagen original (784 dimensiones).
Por lo que es posible "generar" nuevas imágenes eligiendo puntos aleatoriamente del plano 2d, y pedirle al modelo PCA aprendido que invierta la transformación para obtener las "teóricas" imágenes que habrían sido proyectadas a esos puntos del espacio proyectado.
# Mínimos y máximos
c1_min = min(X_pca[:,0])
c1_max = max(X_pca[:,0])
c2_min = min(X_pca[:,1])
c2_max = max(X_pca[:,1])
# Secuencias
c1_sec = np.linspace(start=c1_min, stop=c1_max, num=10)
c2_sec = np.linspace(start=c2_min, stop=c2_max, num=10)
Con las dos secuencias de 10 (una por cada dimensión del plano de proyección) valores es posible combinar los puntos de ambas secuencias para generar 100 combinaciones (puntos 2d) que teselan el plano sobre el que PCA ha proyectado las muestras.
# Método que devuelve las secuencias ordendas en el espacio proyectado
def sorted_sec(c1_sec, c2_sec):
secuencia = []
c2_sec_aux = c2_sec
while len(c2_sec_aux) > 0:
maximo = max(c2_sec_aux)
pos_max = np.argmax(c2_sec_aux)
for c1 in c1_sec:
secuencia.append([c1,maximo])
c2_sec_aux = np.delete(c2_sec_aux, pos_max)
return np.array(secuencia)
secuencia = sorted_sec(c1_sec, c2_sec)
fig, axis = plt.subplots(10, 10, figsize=(10, 10))
row = 0
col = 0
for i in range(len(secuencia)):
if i % len(c1_sec) == 0 and i > len(c1_sec) - 1: # Si empieza una nueva fila y no es la primera
col = 0
row += 1
number = modelo_pca.inverse_transform(secuencia[i])
axis[row,col].imshow(number.reshape(28, 28), cmap='gray_r')
axis[row,col].axis('off')
col += 1
plt.tight_layout()
Sí que se puede interpretar las imágenes reconstruidas, ya que éstas se corresponden con el espacio proyectado anteriormente, pero no son del todo claras dichas representaciones.
Por ejemplo, en la parte superior vemos que el número que más se intuye es el 9 y a medida que nos desplazamos a hacia la derecha ese nueve se va transformado hasta llegar a una especie de 0 sin ser el 0. Si partimos de la zona media de la representación vemos que el número que más se intuye es el 1 y éste se va transformando por 3/5/8 hasta llegar al 0. Y si nos fijamos en la zona baja de la representación vemos como una especie de 1 sin llegar a ser uno, se va transformando por 2 y 3 hasta llegar a una especie de 2 y 0.
Es decir, vemos que sí se generan números y transiciones creíbles según el espacio proyectado que se ha graficado al principio del ejercicio, pero una cosa es creíble y otra entendible, ya que de lo segundo se requiere de una gran abstracción. Esto sucede porque hemos generado diferentes puntos dentro del espacio proyectado, para ser más exactos 100, los cuales de forma equidistante nos permiten a partir de PCA obtener su representación, lógicamente estos puntos no son del dataset original, por lo que pueden quedar bien representados o ser una mera transición entre números.
En resumen, vemos que sí que se consigue obtener números o transiciones no muy bien definidas entre números a partir de puntos equidistantes en el plano gracias a poder calcular la transformación inversa.
No, no se podría haber conseguido lo mismo con t-SNE, esto se debe a que nos permitiría realizar la primera parte del ejercicio pero no la segunda, es decir, t-SNE nos permitiría visualizar las diferentes clases en el espacio proyectado, pero no nos permitiría calcular la transformación inversa, por lo tanto no conseguiríamos que a partir de unos puntos equidistantes del plano se pudieran obtener los números o las transiciones entre los mismos.
# Creación del UMAP
modelo_umap = umap.UMAP(n_components=2, random_state=0)
X_umap = modelo_umap.fit_transform(X)
# Visualización
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(X_umap[:, 0], X_umap[:, 1], s=5, c=y_numerico, cmap='Spectral')
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('Análisis de UMAP', fontsize=14)
plt.xlabel('Primer componente principal', fontsize=14)
plt.ylabel('Segundo componente principal', fontsize=14)
Text(0, 0.5, 'Segundo componente principal')
En este caso UMAP a nivel visual se comporta mejor que PCA, ya que sí que quedan visiblemente separadas cada una de las clases.
Vemos que podemos seguir identificando patrones, ya que hay números que utilizan unos píxeles similares como el grupo formado por los números 7, 9 y 4. Otro grupo formado por los números 8, 5 y 3, y el resto de números que son grupos independientes, aunque por ejemplo el número 6 está más cerca del 0 que del 1.
La principal diferencia entre PCA y UMAP es qu el primero es un algoritmo lineal y el segundo no, es decir, PCA está basado en la factorización de matrices mientras que UMAP está basado en gráfos, es por ello que UMAP consigue identificar mejor a nivel visual las clases.
# Mínimos y máximos
c1_min = min(X_umap[:,0])
c1_max = max(X_umap[:,0])
c2_min = min(X_umap[:,1])
c2_max = max(X_umap[:,1])
# Secuencias
c1_sec = np.linspace(start=c1_min, stop=c1_max, num=10)
c2_sec = np.linspace(start=c2_min, stop=c2_max, num=10)
# Secuencia ordenada para respetar el espacio proyectado
secuencia = sorted_sec(c1_sec, c2_sec)
# Números de la representación
representaciones = modelo_umap.inverse_transform(secuencia)
# Representación
fig, axis = plt.subplots(10, 10, figsize=(10, 10))
row = 0
col = 0
for i in range(len(representaciones)):
if i % len(c1_sec) == 0 and i > len(c1_sec) - 1: # Si empieza una nueva fila y no es la primera
col = 0
row += 1
number = representaciones[i]
axis[row,col].imshow(number.reshape(28, 28), cmap='gray_r')
axis[row,col].axis('off')
col += 1
plt.tight_layout()
En este caso vemos que las imágenes reconstruidas son perfectamente entendibles, más que lo que sucedía con PCA, es decir, nos genera principalmente números aunque también hay alguna que otra transición.
EL porqué es sencillo, UMAP trata de "clasificar" los datos además de reducir la dimensionalidad, con ello me refiero a que busca que los grupos de datos estén lo más diferenciados posible, mientras que PCA busca proyectar cada punto en el hiperplano buscando la máxima variabilidad independientemente de si los grupos están más cerca o no. Esto hace que cuando hagamos la transoformación inversa UMAP funcione mejor, ya que hay una mayor separación/diferencia entre los grupos, mientras que en PCA creaba un grupo muy grande en el que estaban todos los datos.
En este ejercicio se busca automatizar la localización de lugares turísticos a través de los metadatos de las fotografías de flickr.
Para ello se provee junto a la PEC el dataset: barcelona.csv. Ya que se pide encontrar los puntos de mayor interés turístico de esta ciudad.
Opcional: si quieres hacerlo para otra región
Pero si quieres hacerlo para otra parte del mundo, puedes descargarte el dataset completo aquí y descomprime para extraer el CSV.
Para seleccionar las coordenadas de la zona de interés puedes usar la opción Export manual de OpenStreetMaps.
Por último, para filtrar los datos que se corresponden a la zona deseada puedes usar el programa AWK mediante la siguiente línea:
awk -F"," 'NR == 1 {print $5","$6} (NR > 1 && $5 > 41.3560 && $5 < 41.4267 && $6 > 2.1300 && $6 < 2.2319) {print $5","$6}' photo_metadata.csv
$5 hace referencia a la latitud, y $6 a la longitud. Sustituye los valores mínimo y máximo para obtener los datos de localización referentes a tu área de interés.
geo_df = pd.read_csv('./data/barcelona.csv', header=0)
geo_df.sample(5)
| latitude | longitude | |
|---|---|---|
| 12408 | 41.404778 | 2.174541 |
| 2740 | 41.391833 | 2.165167 |
| 1701 | 41.368942 | 2.147912 |
| 11082 | 41.403145 | 2.212886 |
| 4379 | 41.386500 | 2.155999 |
def show_barcelona(df):
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(df.latitude, df.longitude, s=10, alpha=.25)
plt.title("Lugares turísticos de Barcelona")
plt.xlabel("Latitud")
plt.ylabel("Longitud")
show_barcelona(geo_df)
El algoritmo más adecuado para aplicar clustering sobre el problema es DBSCAN, es decir, un algoritmo de densidad.
Es verdad que se podría hacer uso de otros algoritmos como por ejemplo KMEANS, pero éste se basa en el cálculo de centroides y asigna cada punto al centroide más cercano, esto puede ser un inconveniente porque dado este problema podría considerar varios clústers como uno mismo. Además, tanto KMEANS como los algortimos jerárquicos necesitan pasarle por parámetro al algoritmo el número de clústers, cosa que en este ejercicio se desconoce, y lo que se busca es conocer eso número de clústers, es decir, cuántos puntos turísticos hay en Barcelona. También ha que destacar que DBSCAN es el único en detectar outliers, siendo éstos los puntos menos turísticos.
En resumen, el algoritmo DBSCAN es el que mejor funciona según nuestro problema ya que se basa en la estrucutra/densidad de los datos, de esta forma podemos conseguir qué clústers hay en nuestro problema y cómo son.
geo_df_muestra = geo_df.sample(frac=.3, random_state=0)
show_barcelona(geo_df_muestra)
# Estmiación de epsilon
elbow_dbscan(geo_df_muestra, 20)
# Modelo
modelo_dbscan_barcelona = cluster.DBSCAN(eps=0.0025, min_samples=20)
modelo_dbscan_barcelona.fit(geo_df_muestra)
# Visualizamos
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(geo_df_muestra.latitude, geo_df_muestra.longitude, s=10, alpha=.25, c=modelo_dbscan_barcelona.labels_, cmap="tab20")
plt.title("Lugares turísticos de Barcelona")
plt.xlabel("Latitud")
plt.ylabel("Longitud")
Text(0, 0.5, 'Longitud')
# Limpieza, si el clúster es -1 entonces es un outiler
geo_df_muestra["Target"] = modelo_dbscan_barcelona.labels_
geo_df_clean = geo_df_muestra[geo_df_muestra["Target"] != -1]
# Visualizamos
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(geo_df_clean.latitude,
geo_df_clean.longitude,
s=10,
alpha=.25,
c=geo_df_clean.Target,
cmap="tab20")
plt.title("Lugares turísticos de Barcelona")
plt.xlabel("Latitud")
plt.ylabel("Longitud")
Text(0, 0.5, 'Longitud')
Antes de hacer el análisis vamos a calcular el punto medio de cada clúster, el cual nos va a servir tanto para este ejercicio como para el siguiente:
# Método que devuelve la media de la longitud y latitud para cada clúster
def get_meanCluster(df):
dicMeanCluster = {}
for i in np.unique(df.Target):
latitud = np.mean(geo_df_clean[geo_df_clean.Target == i].latitude)
longitud = np.mean(geo_df_clean[geo_df_clean.Target == i].longitude)
dicMeanCluster[i] = [latitud, longitud]
return dicMeanCluster
# Método que imprime la media de la longitud y latitud de cada clúster
def print_MeanCluster(dic):
print("\t\t\tLatitud\t\tLongitud")
for k,v in dic.items():
print("Clúster " +
str(k) +
"\t->\t" +
str(round(v[0], 7)) +
"\t" +
str(round(v[1], 7)))
mean_clusters = get_meanCluster(geo_df_clean)
print_MeanCluster(mean_clusters)
Latitud Longitud Clúster 0 -> 41.3851449 2.1749851 Clúster 1 -> 41.3713111 2.150776 Clúster 2 -> 41.3970314 2.1904491 Clúster 3 -> 41.4132301 2.1991785 Clúster 4 -> 41.4015635 2.1682176 Clúster 5 -> 41.4062836 2.2015134 Clúster 6 -> 41.3714757 2.1402785 Clúster 7 -> 41.4036251 2.1900497 Clúster 8 -> 41.4107145 2.2215941 Clúster 9 -> 41.4038759 2.1743124 Clúster 10 -> 41.4144028 2.1526811 Clúster 11 -> 41.3645843 2.1642466 Clúster 12 -> 41.3753117 2.1570397 Clúster 13 -> 41.3805483 2.1546133 Clúster 14 -> 41.3970029 2.1946244 Clúster 15 -> 41.3984371 2.2124015 Clúster 16 -> 41.3983408 2.1824438 Clúster 17 -> 41.3783561 2.1613572 Clúster 18 -> 41.3688051 2.1608188 Clúster 19 -> 41.4116274 2.1743459
Para interpretar cada clúster se ha calculado el punto medio de cada uno y se ha buscado dicha latitud y longitud en Google (por desconocimiento de la ciudad de Barcelona). Es verdad que para muchos de los puntos calculados sí que aportan sitios significativos, pero para otros no tanto, seguramente haya más que puntos zonas/barrios de interés.
A partir de la anterior visualización y la búsqueda en Google hemos obtenido esta información:
Cabe destacar que la visualización es un poco costosa de entender, esto se debe a que la parte más alta hace referencia al mar, esta parte se puede intuir así, pero si al "norte" está el mar mirando Google Montjuic tendría que estar a la derecha, y vemos que esto no es así, es decir, Montjuic está a la izquierda. Esto nos da a entender que la visualización no es la realidad sino que es un modo espejo de la realidad.
import smopy
mapa = smopy.Map((41.36, 2.13, 41.42, 2.23), z=12)
lista_x = []
lista_y = []
for k, v in mean_clusters.items():
x, y = mapa.to_pixels(v[0], v[1])
lista_x.append(x)
lista_y.append(y)
ax = mapa.show_mpl(figsize=(15,15))
ax.scatter(lista_x, lista_y, s=350, c="purple")
<matplotlib.collections.PathCollection at 0x25e21dc97c0>